subroutine GL(step,energy_out)
use var
use type
implicit none
integer :: iout, tt
integer, intent(inout) :: step
real(dp) :: delta, delta2, dmc_cur, temp_cur
real(dp) :: energy(5), E1, E2, DE, ac, dt_cur, max_ac, min_E
real(dp), intent(out) :: energy_out(5)
logical :: flag = .true.

allocate(q_new(6,tnodes))
iout = 300 + 1

delta = (dt(2)-dt(1))/dt(3)
dt_cur = dt(1)

if (step/=0) then
   call free_energy(energy)
   write(*,'(I8,2X,4(f15.6),2(f20.6))')step,energy,Escale*energy(5)/kbT
   dt_cur = dt(1) + delta*step
   if (step>dt(3)) dt_cur = dt(2)
   call  ginzburg_landau(dt_cur)
   step = step + 1
endif
E1 = 0   ! Initial Energy
ac = 1   ! accuracy
tt = 1   ! parameter for trj_E
flag = .true.
Do while (flag)
   if (mod(step,sample)==0) then
      call free_energy(energy)
      ! calculate accuracy; energy between two step
      E2 = energy(5)
      DE = (E1-E2)
      ac = abs(DE/E1)*100
      E1 = E2

      write(iout,'(I8,2X,4(f15.6),2(f20.6))')step,energy,Escale*energy(5)/kbT
         open (22, file = 'timestep')
         write(22,'(I8,2X,4(f15.6),3(f20.6))')step,energy(1:5),Escale*energy(5)/kbT,Escale/kbT
         write(*,'(I8,2X,4(f15.6),3(f20.6))')step,energy(1:5),Escale*energy(5)/kbT,Escale/kbT
         close(22)

      trj_E(:,tt) = energy
      tt = tt+1
      if (ac < accuracy(2)) flag = .false.

      dt_cur = dt(1) + delta*step
      if (step>dt(3)) dt_cur = dt(2)
   endif
   if (mod(step,traj(1))==0) then
      write(*,'(I8,2X,4(f15.6),2(f20.6))')step,energy(1:5),Escale*energy(5)/kbT
      write(222,'(I8,2X,4(f15.6),2(f20.6))')step,energy(1:5),Escale*energy(5)/kbT
      if (traj(2)==1) call wtraj(1,step,step)
      tt = 1
   endif
   
   call  ginzburg_landau(dt_cur)
   
   step = step + 1
enddo
write(222,'(I8,2X,4(f15.6),2(f20.6))')step,energy(1:5),Escale*energy(5)/kbT
energy_out = energy
deallocate(q_new)
endsubroutine GL
